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Abstract 

The quasi-steady state approximation and time-scale separation are commonly applied 
methods to simplify models of biochemical reaction networks based on ordinary differential 
equations (ODEs). The concentrations of the “fast” species are assumed effectively to be 
at steady state with respect to the “slow” species. Under this assumption the steady state 
equations can be used to eliminate the “fast” variables and a new ODE system with only 
the slow species can be obtained. 

We interpret a reduced system obtained by time-scale separation as the ODE system 
arising from a unique reaction network, by identification of a set of reactions and the corre¬ 
sponding rate functions. The procedure is graphically based and can easily be worked out 
by hand for small networks. For larger networks, we provide a pseudo-algorithm. We study 
properties of the reduced network, its kinetics and conservation laws, and show that the 
kinetics of the reduced network fulfil realistic assumptions, provided the original network 
does. We illustrate our results using biological examples such as substrate mechanisms, 
post-translational modification systems and networks with intermediates (transient) steps. 

Key words: Reduced network, quasi-steady-state, species graph, noninteracting, dy¬ 
namical system, positivity 
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1 Introduction 

Biochemical reaction networks often involve many biochemical species that interact through 
many reactions. The mathematical models that are used to describe such networks can be 
quite complex and analytically intractable, both in terms of variables (species concentrations) 
as well as (unknown) parameters such as reaction rate constants. It is therefore commonplace 
to reduce the dimensions of the models by various means. 

One way to simplify a model is by time-scale separation mu- To apply time-scale 
separation, the species in the system are divided into fast and slow species. Fast species are 
assumed to be at equilibrium (a so-called quasi-equilibrium) even though the system as a whole 
has not reached a state of equilibrium. For example, if degradation of a species proceeds at 
high rate, then it is short lived and becomes a ‘fast’ species. It is thus reasonable to assume 
that the species is in a state of (quasi-)equilibrium. Being at equilibrium, the fast species 
might be eliminated from the system, resulting in a simpler, reduced system with only the 
slow species. Tikhonov’s theorem mm might be used to conclude that the original and the 
reduced system have similar dynamics (over a compact time interval). 
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Our aim is to interpret the reduced system as an ODE system arising from a reaction 
network with only the slow species. We follow some ideas in The starting point is a 

system of ordinary differential equations (ODEs) arising from a reaction network with a given 
kinetics. We assume that a set of species (‘the fast species’) is given and that these species 
effectively are at steady state. If the species do not interact with each other (i.e. they are not 
on the same side of a reaction), then a relationship between the concentrations of the fast and 
the slow species can be derived under certain conditions [3]. Using this relationship, a reduced 
ODE system for the slow species in which the fast species have been eliminated is obtained. 

The next step is to interpret the ODE system for the slow species as the ODE system 
associated with a reduced reaction network. This step could be carried out in various ways. 
We identify sequences of reactions in the original network, by means of a specific graph (see 
example below), such that the net production of the fast species is zero. By contracting the 
sequence into a single reaction and removing the fast species, a reduced reaction is obtained. A 
non-trivial issue is to determine the reaction rate of the reduced reaction. Eor small networks, 
the procedure can easily be carried out by hand. 

Subsequently, we establish that the reduced reaction network has some desirable basic 
properties. Eor example, if the reactions of the original reaction network cannot take place in 
the absence of their reactant species, then this also holds for the reduced reaction network. 

Graphical means to re-interpret a reaction network under quasi-stationarity are not new. 
The King-Altman procedure is a systematic way to eliminate enzymes and intermediate 
species, in an enzyme catalysed system with mass-action kinetics [12j . Later, Wong and 
Hanes [23] gave a systematic method to find the production rates of the slow species, avoiding 
much of the algebraic manipulations in [12]; see also [9]. These approaches can be seen as 
instances of a general linear (graphical) elimination procedure [3] 0] . 

None of the mentioned approaches specifically identifies a reduced reaction network, except 
for [ 1 ], where a reduced reaction network is obtained after elimination of intermediate species. 
A related approach by Horiuti, Temkin and co-workers, is based on a graphical procedure to 
identify, not one but many, reduced reaction networks [mn]. The procedure is similar in 
spirit to our method and we will contrast the two approaches in Subsection 17.11 In [7] an 
algorithmic method to obtain possible reaction networks corresponding to a rational ODE 
system is described. In m a graphical method for the elimination of complexes (not species) 
is presented. 

To illustrate the results and methods of the paper, consider a ping-pong bi-bi mechanism 
[ 6 ], described by the reaction network 

E + Si^Yi^E* +Pi E*+ S2 ^Y2^E + P2, 

k2 ki fcg kg, 

where E, E* are two forms of an enzyme. Si, S 2 are substrates, Yf, Y 2 intermediates and Pi, 
P 2 products. The molar concentrations of the species are denoted as: xi = [P], X 2 = [E*], 
X 3 = [Pi], xi = [P 2 ], X 3 = [Pi], xq = [P 2 ], X 7 = [Yi] and xg = [Y 2 ]. With this notation and 
assuming mass-action kinetics, the evolution of the species concentrations is described by the 
ODE system (cf. Equation ([ 6 |)): 


xi = k2Xf + kjxs - {kiX3 ksXQ)xi 
X3 = k2Xj - kiXiX3 
X5 = /C3X7 - ^4X2X5 
X7 = A:iXiX3 -I- k 4 ,X 2 X 3 - {k2 -I- k3)x7 


X2 = /C 3 X 7 -I- kexs - {kiX3 + k3Xi)x2 
X4 = kexs — ^5X2X4 
xe = kjxg — kgxiXQ 
xg = ^ 5 X 2 x 4 -I- kgxixe - (ke + k7)xg 
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The main interest is on the conversion of substrates into products [2j- It is thus sensible 
to eliminate the set of species U = {E,E*,Yi,Y 2 } to obtain a reduced reaction network on 
the species 5i, S 2 , -Pi, P 2 alone. Assuming that P, E*,Yi, I 2 are at steady state and using the 
conservation law for the total amount of enzyme, T = xi + X 2 + xj + we obtain 

xi = q{x){{kQ + kj)k2k4,x^ + {k 2 + A:3)fc5^7a^4) 

X2 = q{x){{kQ + k7)kik3X3 + (^2 + k3)kGk8Xe) (1) 

X 7 = q{x){kik5k7X3X4 + (fce + kj)kik4X3X3 + kikekgXGXG) 

X8 = q{x){kik3k5X4X3 + k2k4^ksX5X6 + {k2 + k3)k5k8X4Xe), 


where 


q{x) = T j [kik3{k3 + k7)x3Xi + (feg + k7){kik^X3X3 + kiksxs + k2kiX3) 
+ kik8{k2 + kQ)x5X6 + {k2 + A:3)(/c5/c8X4X6 + k3k7Xi + keksXe)^. 


Conditions that guarantee existence and positivity of functions expressing the concentrations 
of the species in a set U in terms of the remaining species, were given in [3]. These are reviewed 
in Section [3l The form of the expressions is the content of Theorem [11] 

After substitution of ([I|) into the ODE system, we obtain the production rates of Pi, P 2 
and the consumption rates of Pi, 82 - 


X5 = xg = -X3 = -X 4 = q{x){-k 2 kikGk 8 X 3 XG + kik3k3k7X4X3). (2) 


This ODE system might be interpreted as arising from a reaction network with one reversible 
reaction with the following rate functions: 

Pi + P 2 Pi + P 2 , Ki = q{x)kik3k5k7X3X4, K 2 = q{x)k2k4kGk8X5X6. (3) 

K 2 


In this particular case, it seems straightforward to identify reasonable reactions that explain 
(ED. However, in general this might not be so. A procedure to systematically find a reduced 
reaction network is presented in Section 01 At the core of the procedure is a graph that relates 
the species in U (Dehnition [7D. The graph of the example is 



Each edge in the graph corresponds to a reaction in the original network. Eor example, the 
edge E — Y 2 corresponds to the reaction P + P 2 12 - 

Our main result is Theorem 1161 that establishes the reduced reaction network. Essentially, 
it has two types of reactions: The reactions of the original network that do not involve species 
in U, and reactions found by considering certain cycles of the graph. The reactant (resp. 
product) of such a reduced reaction is the sum of the reactants (resp. products) of the original 
reactions defining the cycle, after removing the species in U. Eor example, the above graph has 
two cycles that give rise to two reduced reactions. The clockwise cycle involves the reactions 
corresponding to the labels kiX 3 , /cs, ^ 5 X 4 and fey: 

E + Si^Yi Yi%E* + Pi E* + S 2 ^Y 2 Y 2 %E + P 2 . 
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Adding the reactants and products together and removing E, E* , Yi, and Y2, we obtain the 
reaction Si + S 2 ^ Pi + P 2 - Similarly, the anti-clockwise cycle gives the reaction Pi -|- P 2 ^ 
Si + S 2 (consistent with ([3D). The rate function of each reaction is determined from the edge 
labels of the graph. We describe an algorithm to find the reduced reaction network and the 
rate functions after Theorem 1161 

Note that the reduced rate functions in ([3D have ‘mass-action form’ in the sense that the 
reactions can only occur in the presence of the reactant species. This holds generally for the 
reduced reaction network (cf. Section [5D . 

The outline of the paper is as follows. In Section [2] we introduce background material. 
In Section [3| we recall some results concerning elimination of variables in [3] . The reduced 
reaction network is derived in Sectional Properties of the reduced reaction network in relation 
to the kinetics (including mass-action) and conservation laws are presented in Section [5j This 
concludes the core of the paper. Sections [Gj [7| and [8| discuss respectively iterative elimination, 
comparison with previous results (the approach in [T7| and networks with intermediates) and 
post-translational modification networks. Finally, Section [9| contains proofs. 

2 Preliminaries 

In this section we introduce necessary concepts from graph theory and reaction network theory. 

We let M>o and M>o denote the sets of nonnegative and positive real numbers respectively, 
and define M>q and accordingly. For x, y G M”, x ■ y denotes the scalar product associated 
with the Euclidean norm. Further, {vi, ... ,Vr) denotes the vector subspace generated by 

Vi, ... ,Vr € M"’ . 

Graphs, multidigraphs and spanning trees. Let Q = be a directed graph 

(digraph) with node set M and edge set £. By abuse of notation, we write e G ^ whenever 
e G A spanning tree r is a directed subgraph of Q with node set J\f and such that the 
underlying undirected graph is connected and acyclic. A spanning tree is rooted at the node 
if is the only node with no outgoing edges. 

The graph G is strongly connected if there is a directed path from A^i to N 2 for any 
pair of nodes A^i, N 2 . Any directed path from A^i to N 2 in a strongly connected graph can 
be extended to a spanning tree rooted at N 2 . A cycle is a closed directed path Ni^ —>■ 

■ ■ ■ ^ Ni^ —>■ with no repeated nodes apart from the initial and terminal nodes. By 

definition all cycles are directed. 

If TT i? is a labeling of Q with values in some ring R, then any subgraph % oi G 

inherits a labeling from G- We extend the function vr to the set of subgraphs of G by defining 

^('^) = n 

eeH 

A multidigraph ^ is a pair of finite sets {N',£) equipped with two functions: 

s:£^N t:£^M. (4) 

The elements of N are called nodes, the elements of £ are called edges, and the functions 
s,t are the source and target function, respectively. The function s assigns to each edge the 
source node of the edge and the function t assigns to each edge the target node of the edge. 
In a multidigraph both self-edges (edges e with t{e) = s(e)) and parallel edges (edges 61,62 
with t(ei) = t(e 2 ) and s(6i) = 5 ( 62 )) are possible. 
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Spanning trees, cycles and labels for a multidigraph are defined analogous to those of a 
digraph. Note that the unique spanning tree of a multidigraph with one node and one self-edge 
is the node itself. 

We associate a digraph Q with a multidigraph Q by removing self-edges and collapsing 
parallel edges into one edge. That is, only one of the parallel edges between two nodes is kept. 
If Q is labeled, then so is Q and the label of an edge in Q is the sum of the labels of the parallel 
edges in Q with the same source and target. If a node in G is only connected to itself then it 
is not included in A formal definition of Q is given in Subsection 19.11 

Reaction networks. A reaction network on a finite set 5 is a multidigraph (C, TZ) 
where 

(i) S = {5i, ..., S'n} is called the species set. It is equipped with an order that provides a 
canonical isomorphism = M”. 

(ii) C C M>g is called the set of complexes. 

(hi) TZ = {ri, ..., ri} is called the set of reactions. 

The source function in (Uj) is denoted as y: TZ ^ C and assigns to each reaction r* its reactant 
yr^ ■ The target function in (jU is denoted as y': TZ ^ C and assigns to each reaction r* its 
product ?/(.,. 

Elements of are identified with linear combinations of species. Hence, under the iso¬ 
morphism = M”" in (i), we write a complex as a linear combination of the species. For 
instance, the complex (1,0,1) G C, which corresponds to the element in assigning 1 to Si 
and S 3 and 0 to S'2, is written as 5 i -|- ^s. 

It is assumed that yr 7 ^ for all reactions r ^TZ, and that every complex t] £ C is either 
the reactant or the product of some reaction r £TZ. That is, the multidigraph (C, TZ) contains 
no self-edges and no isolated nodes. 

Given a complex rj £ C G M>q we call ry the stoichiometric coefficient of Si in rj. We 
say that y involves Si ii rji 7 ^ 0 , and that a reaction r involves Si if Si is involved in the 
reactant or product of r. We say that a pair of species Si,Sj £ S, i ^ j, interact if they are 
both involved in the same reactant or product of a reaction. Equivalently, if r]i,rij 7 ^ 0, for 
some complex y. 

It is not common to allow for multiple reactions between the same reactant and product. 
It is however mathematically convenient for our purposes, as will be clear in Section 01 Note 
that, if ^ is a reaction network in our sense, then ^ is a reaction network in the standard 
sense. 

We often give a reaction network by listing its reactions. The set of complexes C and the 
set of species S is easily found from the reactions. 

Example 1. The following 

Si + 54 ^ 55 —7 52 + 54 —7 S3 + S4 

is a reaction network on the set of species S = {Si, S2, S3, S4, S5} with set of complexes 
C = {Si+ S4,S5,S2 +54,53 + 54}. 
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Dynamical systems. A kinetics for a reaction network (C, TZ) is a function 

KlO—>-]R>g X !->■ (Kri(x), . . . , Kr^(x)), 

where ]R”q C 0 C ]R>g, such that k(M”q) C M^q. The component i^nix) of k{x) is called the 
rate function of the reaction r* G TZ. 

Mass-action kinetics is defined by the following rate functions: 

n 

Kr{x) = krX^'" = kr r ^ TZ, 

i=l 

where kr > 0 is called the reaction rate constant of the reaction r. By convention, 0*^ = 1. 

The kinetics k provides a labeling of the reaction network. Since the set of reactions is 
ordered, we often denote the rate functions as Ki{x) instead of i^nix)- In the examples, we 
label reactions with their rate functions. In the particular case of mass-action kinetics, we 
simply use the reaction rate constants ki as labels. 

We denote by Xi the concentration of species Si and let x = (xi,..., Xn) be the vector of 
concentrations. In specific examples, the species are denoted by letters such as E, P and the 
concentrations are denoted by xe,xp, respectively. 

For a reaction network {C,TZ) on S and a kinetics k{x), we let 

9i{x) ='^ Kr{x){yr - yr)i, i = X e fi. (5) 

The evolution of the species concentration in time is modelled by the following system of 
DDEs: 

X = g{x), X G n, ( 6 ) 

where x = (xi,... ,Xn), g{x) = (< 71 ( 3 :),... , 5 „(x)), and where the derivative is with respect to 
time. Explicit reference to time is omitted. 

The steady states of the system ([ 6 |) are the solutions to the system 

g{x) =0, X G n. 

Remark 2. Consider a reaction network Q with a kinetics k and the associated reaction 
network Q with the induced kinetics k, (the induced labeling of G). Then the ODE system 
associated with {G,k) agrees with the ODE system associated with {G,k). 

Remark 3. In typical models of biochemical reaction systems, the rate function of a reaction 
vanishes whenever the concentration of one of the reactant species is zero. In particular, this 
guarantees invariance of the non-negative orthant under ([ 6 ]). However, we do not need this 
assumption for our results to hold. In Section 15.11 we discuss some results that follow from 
making this assumption. 

Conservation laws. The stoichiometric subspace of a network (C, TZ) is the vector 
subspace of M” given by 

S={y'r-yr \ ren)c M”. 

If a; = (cji,... ,Un) G S'"*", then it follows from ([5]) and ([ 6 ]) that uj ■ x = 0. Thus, for any 
trajectory there is a constant T G M such that 

n 

T = u ■ x = uJiXi. 

i=l 
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This equation is called the conservation law with total amount T G M, corresponding to 
u G S-^. We also say that w • x is conserved. A set of conservation laws 

{Ti = • X, ... , Ti = iJ ■ x} 

is minimal if form a basis of S-^. Then, the trajectory with initial concentration 

xo is confined to the linear space with equations 

Tj = cj* • X, with Tj = w* • Xo, for i = 1,..., /. 

Example 4. Consider Example [T] and let k be a kinetics: 

s. + S4 ^ S, + S4 ^ S3 + S4. 

K2(x) 

The corresponding ODE system is: 

Xi = - Kl(x) + «: 2 (a^), X2 =K3{x) - K4{x), X 3 =K 4 (x), 

X4 = - Ki(x) + /t2(x) + «:3(x), X5 =Ki{x) - K2 {x) - K3(x), 

and a minimal set of conservation laws consists of 

Xi + X 2 + X 3 + X 5 = Ti, X4 + X 5 =r 2 . (7) 

Since concentrations are nonnegative, Ti,T 2 > 0. 

3 Elimination of variables 

In this section we introduce some results from [3] about linear elimination of variables. The 
goal in [3] is to express the concentrations of some of the species at steady state in terms of 
the concentrations of the other species. This is done for sets of species, say U, that form a 
noninteracting set (Definition [5]) , assuming the kinetics is linear with respect to the concentra¬ 
tions of the species in U (Definition [B]). The elimination makes use of a special multidigraph, 
called Qk , which is defined in Definition [7] (see also Theorem ITT]) . This multidigraph is also 
used to define the reduced reaction network in Section [H 
Let a reaction network (C, TV) on a set S be given. 

Definition 5. A subset U C S is noninteracting if it contains no pair of interacting species, 
and the stoichiometric coefficients of the species in U in all complexes are either 0 or 1 . 

Let 7/ C 5 be a noninteracting subset of species. Eor simplicity we let U = {Ui ,..., Um} 
and = S \ U = {S*!,..., Sp} (with p = n — m), such that S = UU. We order S as 
S = {5i,..., Sp, Ui ,..., Um}- Erom now on, we let x* be the concentration of Si G and 
X = (xi,... ,Xp). Similarly, we let Ui be the concentration oi U ^ U and u = (tti,... ,Um)- 
Hence, the rate functions are functions of (x,u): Kr{x,u). 

We let 

C:R^^RP, ( 8 ) 

be the projections onto the first p coordinates and last m coordinates of M”, respectively. 
Eurther, we let TZu be the set of reactions that involve species in U in the reactant and/or in 
the product: 

TZu = {r en I p{yr) / 0 or p{y}) 7 ^ 0}. 
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Any reaction in TZu involves at most one species in lA in the reactant and at most one in the 
product. Hence, the vectors p{yr) and p{yr) have at most one nonzero component (equal to 
one) for any r € TZu- 

We next impose some regularity conditions on the rate functions, namely, that the rate 
functions are linear in the concentrations of the species in Zi in a specific way. We assume H 
takes the form H x M>g, where M^g C H C M^g (making an abuse of notation). 

Definition 6. Let k: H x M>q ^ M>g be a kinetics with M^g C H C M>g. The kinetics k is 
^//-linear if, for each r € TZu, there exists a function H —>• ]R>o such that Vr{x) > 0 for all 
a: € H and 

, . / UiVr{x) iir gTZu and p{yr)i = 1 

''’'<^•“> = 1 vAx) it . e and pfe,) = 0. 

A Z//-linear kinetics can be interpreted as being mass-action with respect to the species in 

U. 


Definition 7. Let Z// C 5 be a set of noninteracting species and k a Z^Z-linear kinetics. We 
define the labeled multidigraph Qu = {Nu-,£u) by 


and 




U if p{yr) / 0 and p(y(.) / 0 for all r E TZu 

U\J {*} otherwise 


£u = {Ui Uj I r E TZu with p{yr)i / 0 and p{yr)j / 0} U 

{Ui * \ r eTZu with p{yr)i 7 ^ 0 and p(y(,) = 0} U 

{* Ui \ r eTZu with p{yl.)i / 0 and p{yr) = 0 }, 

where the last two subsets in the definition of £u are empty when Mu = U. 


Observe that Qy might contain parallel edges between any pair of nodes and self-edges for 
nodes other than *. Each edge e in the multidigraph Qy corresponds to a reaction in TZy, 
r(e). Moreover, since U is noninteracting, the multidigraph Qy has exactly one edge e(r) for 
each reaction r in TZy. This gives rise to two bijective functions 


£u 



(9) 


such that r o e = idT^^^ and e o r = idf^^. In the examples, the functions r, e are implicitly 
given by the subindices of the functions Vr{x): the edge label Vi{x) indicates that the edge 
corresponds to the reaction r*. 

The graph Gy (not taking into account the labelling) agrees with the embedded network 
of {C,TZ) given hy U, as defined in |TT] considering * as the zero complex. 

Example 8 . Consider the network 

Ki(x,u) K2 {x,u) K3(x,u) o < tt „ 

- > 02 -^ tVi - > U2 03 + Ui - > U3 -f- 01 

S 3 + (7, S, + t /3 f/3^^S3. 

The set U = {Ui, U 2 , U 3 } is noninteracting. A W-linear kinetics k{x,u) takes the form 
k{x,u) = {Ki{x,u),V 2 ix),UiV 3 {x),UiV 4 {x),U 2 V 5 {x),U 3 V 6 {x)) 
and the corresponding multidigraph Qy is 
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V4{x) 



Example 9. Consider a phosphorylation mechanism for a substrate S such that phosphory¬ 
lation may be catalyzed by two different enzymes, Ei and E 2 , with mass-action kinetics: 


S + Ei^Yi%Sp + Ei S + E2^Y2% Sp + E2. 

k2 k5 

We consider the noninteracting set U = {Ei, E 2 ,Yi,Y 2 }. Since we assumed mass-action ki¬ 
netics, the kinetics is ZY-linear. The multidigraph Qk is 





Define the support of a vector a; = (wi,..., ujn) € M" as 

supp(a;) = {Si I uji ^ 0} Y S. 

A subset T-L YU oi a, noninteracting set is itself a noninteracting set and a Z^Z-linear kinetics 
is 'H-linear. Hence the multidigraph Q-^ is defined. 

Any connected component of Qu agrees with for some T-L YU. If * 0 My, then the 
reactant and the product of any reaction in TZy involve each exactly one species in "H. In this 
case there is a vector oM € S'^ such that 

Q{uM) = Q and p{uM)i = \ ' z = l,...,m, (10) 

that is, the coordinates of uM corresponding to the species in % are all one and only these are 
non-zero. 

According to [3] , any vector in S’-*- with support in ZY is a linear combination of the vectors 
uM for all connected components Qy such that * 0 My. Since these connected components 
are disjoint, the vectors uM are linearly independent. These vectors give rise to conservation 
laws: 

Ui = Ty. (11) 

i I Ui&V. 

In Example [9l there are two connected components without the node *. They give rise to 
two conservation laws, namely, xei+ xy^ = Ti and xe 2 + xy 2 = T 2 . Any other conservation 
law involving only the concentrations of Ei, £' 2 , Ti, T 2 is a linear combination of these two 
conservation laws. The only connected component of Gu in Example [8] contains *. Hence 
there is no conservation law involving only concentrations of species in U. 

Let iV be a node of Gu and Gy be the connected component that contains it, that is 
N G My. Let Q{N) be the set of spanning trees of Gy rooted at N. 

Definition 10. Let (C, TZ) be a reaction network on a set <S and k a kinetics. A set U Y S is 
said to be linearly eliminable if 
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(i) lA is noninteracting. 

(ii) For every connected component of Qu it holds: 0(*) 7 ^ 0 if * € Ny and UArg-^0(iV) 7 ^ 
0 if * 0 J\fy. 

(hi) K is ^-linear. 

Condition (ii) states that Qy admits a spanning tree rooted at *, if * G My, and rooted at 
an arbitrary node, if * 0 My. Any strongly connected component of Qu fulfills this condition. 

Let U Q She a. linearly eliminable set. For each connected component Qy of Qu, let qy{x) 
be the function 


qn{x) = 

Ty 

if * 0 My, where Ty G M>o, 

( 12 ) 

E <ry 

N&n,T&0{N) 




Qn{x) = 

1 

if * G My. 

(13) 

E vr(r)’ 


tS0(*) 

The function qy{x) is a function of x since 7 r(r) depends on x. We omit the explicit reference 
to Ty in qy{x) and to x in 7r(T) for convenience. When the multidigraph Qu is connected, we 
omit the subscript % = U. 

Definition frUTiil guarantees that the denominators are never empty sums. By positivity 
of Vr{x), the function qy{x) is positive for all x However, the denominator of qy{x) 

might vanish at points at the boundary of H. Henceforth, we let H be the domain where qy{x) 
is well defined for all connected components Qy. We have C H C H. 

We will make use of O Propositions 8.4, 8 . 6 ], but reformulate them in the terminology of 
this paper (see also Section [92]) • 

Theorem 11 (Elimination [3|). Let {C,TZ) be a reaction network on a setS, k a kinetics and 
assume U Q S is a linearly eliminable set. Consider the system of equations 

u = 0, and Ty = Ui, Ty > 0, (14) 

for aim QU such that Qy is a connected component of Qu with * ^ My. Equation (fTTl) holds 
for {x,u) € Q X M>q if and only if 

Ui = qy{i){x) ^ 7 r(r), for all i = 1.. .m, 

Tee(Ui) 

where qy(^i-j{x) is given by (fT2]) or (IT^ . and LLii) Q U is such that Ui € ^(i) for Qy(i) a 
connected component of Qu. 

The theorem states that if the species in U are at steady state with respect to the species 
in U‘^, then the concentrations Ui can be given as rational functions in the labels of Qu, for 
fixed total amounts Ty. The labels are independent of the concentrations of the species in U. 
In virtue of the non-negativity of the labels of Qu, the concentrations of the species in U are 
also non-negative for x & Ql. 

For future reference, we define the function (p: H 

<f{x) = {<pi{x),...,<pm{x)), <fi{x) = qy^Pi{x)Y^ 7r(r). (15) 

T&eiUi) 
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If a connected component of Qu is strongly connected, then there is at least one 
spanning tree rooted at each node, and hence (^(x) is not identically zero in any coordinate 
corresponding to a species in %. If is not strongly connected, then at least one node is not 
the root of a spanning tree, and hence at least one component of ^p{x) is identically zero. In 
the latter case, all steady states lay at the boundary of the positive orthant. 

Example 12. Consider the network in Example [T] and the noninteracting set lA = { 54 , 55 }. 
To adapt the notation to the present setting, we let Ui = 54 and U 2 = S^. We assume the 
network is endowed with a Z//-linear kinetics 

k{x,u) = [uiVl{x),U2V2ix),U2V3{x),UlV4^{x)) . 


The multidigraph Qu is 


V4(x) 

Q- -yila;) 


Ui 



V 3 {x) 


The multidigraph Qu is strongly connected and thus the set Li is linearly eliminable. We apply 
Theorem m to conclude that ui = 0 , 'I 62 = 0 and ui + U 2 = T 2 hold if and only if: 


ui = q{x){v3{x) + V2{x)), U2 = q{x)vi{x), q{x) 


T 2 _ 

Vi{x) + V2ix) + V3{x) ' 


Example 13. The set U = {Ui,U 2 ,U 3 } in Example [8] endowed with ZY-linear kinetics is 
linearly eliminable since Qu is strongly connected. By Theorem fTTl at steady state it holds 
that: 


ui = q{x)v2{x)v5{x)v6{x), 

U 3 = qix)v2ix){v3{x) + Vi{x))v5{x), 


U2 = q{x)v2{x)v3{x)ve{x), 

{V3{x) + V4ix))v5{x)vQ{x) ' 


4 Reduced reaction network 

Assume that a reaction network (C, TZ) on a set 5, a kinetics k and a linearly eliminable 
set U Q S are given. If u = (p(x), as given in (fT^ . is substituted into the ODE system 
(x, ii) = g{x, u), given in ([U]), then tt = 0 by construction and an ODE system in the variables 
xi,... ,Xp is obtained: 

X = g{x), X € D C M>g, (16) 

where 

9 {x) = C{9{x,^{x))). 

That is, for i = 1,... ,p, 

h{x) = 9i\xi,...,Xp, qH{i){x) ^ -K{T),...,qu{m){x) ^ vr(T) 

\ re0(;7i) reQiUru) 

In this section we prove that the ODE system (I16p is the ODE system associated with 
a reaction network on the species set {5i,...,5p}. The reactions and the kinetics of this 
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reaction network are graphically obtained from the multidigraph Qu, following the procedure 
given below Theorem 1161 

The intuition is as explained in the introduction. Cycles in the multidigraph correspond 
to sequences of reactions such that the net production of the species in U is zero. However, 
not all cycles give rise to a reduced reaction. In order to give the precise formulation of the 
reduced reaction network (Definition [151 Theorem 1 161) . we will need the following. For a cycle 
cr in let be the connected component that contains the cycle. In particular, the 

nodes of a form a subset of ■H(cr). Let r((T) be the set of subgraphs G of Qu such that a C G 
and such that removing any edge e G a from G creates a spanning tree of G-hG)-: rooted at the 
source s{e) of e. Specihcally, given an edge e € a, dehne 

r((T) = {rUe I a\e G T, r€ 0(s(e))}. (17) 

The set r((T) does not depend on the chosen edge e (Proposition Hi] in Section [931) • We further 
define the function 

n(o-) = TT{j)=Tr{a) vr( 7 \CT), (18) 

7€r((T) 7Gr((7) 

which is the sum of the labels of the graphs in r(iT). 

Remark 14. To hnd r(iT), it is sufficient to consider any edge e G a and hnd all spanning 
trees rooted at s(e), containing cr\e. Therefore, r((T) 7 ^ 0 if and only if the path a \ e can 
be extended to a spanning tree of G'u{a) rooted at s(e). In particular, if is strongly 

connected, then r(iT) 7 ^ 0 . 

If cr is a cycle that contains all nodes of G'H(a)-, then r((T) = {ct} and thus n((T) = 
Indeed, in this case u \ e is the unique spanning tree of GhG) containing cr \ e. 

Let A be the set of cycles a in Gu such that 

- yr{e)) / 0 and r(o-) / 0 . 

eG(J 

The first condition states that the net-production of some species in is non-zero over the 
reactions associated with the edges of the cycle. 

Definition 15. Let {C,TZ) be a reaction network on a set S, k a kinetics and lA G S a 
linearly eliminable set. The rednced reaction network obtained by elimination of lA 
is the reaction network (C,77.) on the species set lA^ with kinetics k defined on D, such that 
7^ = 7^1 U 7 ^ 2 ) where 

^1 = {r: C(yr) C{y'r) I r e 7^ \ Uu) , 

'1^2= ya- ^'^^(y'r(e)) I O' G A 

1 eSc egtr 

C is the set of source and target nodes of TZ, 

Kr{x,(f{x)) if f = C(yr)C(2/0 ^ for r G 7^ \ 7^^./ , 
7-H(<7)(a^)n(<T) ifr = r^G'R2- 

and D C D is the maximal set for which q-uix) is well defined for all connected components 

Gn of Gu- 
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In short, TZi consists of the reactions that do not involve species in U and 7^2 consists of 
new reactions defined by the cycles in A. The procedure outlined below Theorem 1161 clarifies 
how the reactions in 7^2 are obtained. 

We are now ready to present our main result: the interpretation of the ODE system (jl 6 p as 
the ODE system of the reduced reaction network obtained by elimination of lA defined above. 
The proof is given in Subsection 19.31 

Theorem 16 (Reduced reaction network). Let {C,1Z) be a reaction network on a set S, 
K a kinetics and lA ^ S a linearly eliminahle set. Then the ODE system in dm), 

Xi = giix,y:>ix)), X G D C i = l,...,p, 

is the ODE system associated with the reduced reaction network obtained by elimination oflA. 

The fact that we can find a reaction network for the reduced ODE system implies that the 
system does not have negative cross effects in the sense of |16] . 

Procedure. We provide the following procedure to determine the reduced reaction network 
obtained by elimination of a set of species lA, in virtue of Definition [TJ] and Theorem [TUJ 

(1) Check whether lA is linearly eliminable: 

(i) Check that the set lA is noninteracting and that k is W-linear. 

(ii) Compute the multidigraph Qu (as in Definition [7]) and check that each connected 
component Q-n of Qu admits a spanning tree rooted at an arbitrary node, if * is not 
a node of Q-p and rooted at * otherwise. 

If either (i) or (ii) fails, stop. 

( 2 ) ^Compute qn{x) for each connected component Q-^ of Gu and p{x), using (fT 2 ]i . (fT3]l and 

m- 

(3) Compute the reactions in TZi: for each reaction that is not in TZu do 

(i) Add the projection of the reaction by C to TZi. 

(ii) *Compute its rate function by replacing u by p(x) in the original rate function. 

(4) Compute the reactions in 7 ^ 2 : for each cycle a in Qu do 

(i) List the reactions corresponding to the edges in the cycle and add the species in lA^ in 
the reactants and products to obtain the reactant and product of the new candidate 
reaction rv. 

(ii) If y^^ 7 ^ y'~^ proceed. Otherwise go to the next cycle. 

(iii) Compute r((T). If r((T) 7 ^ 0, add r„ to 7^2 and proceed. Otherwise go to the next 
cycle. 

(iv) *The rate function of the new reaction is ( 7 -^( 0 -)( 3 j)n(( 7 ), cf. (fTHjl . 

If we are only interested in the reactions of the reduced reaction network and not the 
kinetics, then the steps marked with * can be ignored. Eurther, in this case, if G-u(a) is 
strongly connected, then r((T) 7 ^ 0 and we do not need to find r(fT) explicitly. 
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Remark 17. Let ri-^,... ,rii € 7^ be reactions in a reaction network (C, TZ) that form a cycle 
in the network. Then, 

l 

Y. ^y'n, - Vn,) = 0 - 

i=i 

Hence, step (|4iH) is not fulfilled and the cycle does not belong to A. This applies in particular 
to any cycle arising from a reversible reaction rji rj 2 . 

Example 18. For illustration, we apply the method to the reaction network in Example [H 
with U = {S' 4 , 55 } = {Ui, U 2 ] and a 77-linear kinetics. 

(1) lA is linearly eliminable: see Example 1121 

(2) q{x) and ^p{x) are given in Example fT^ 

(3) TZii = 'Tl. Therefore 7^i = 0. 

(4) We compute the reactions in 7^2- Qu has three cycles: 

(i) n: Si+JJ{ —^ 

ra: K —> S2+JX 

r.: Si —^ S2 

(ii) = Si ^ 82 = y~y 

(iii) r((Ti) = {(Tl} / 0 (see also Remark [T411 . The reaction 
Si —S 2 belongs to 7^2- 

(iv) Rate function: q{x)vi{x)v'i{x) 



V3(x) 


<T 2 : 


V 2 (x) 


U2 


(i) n: Si+jx —> K 

r2: K —> Si+p{ 

P. Si —^ Si 

(ii) = Si = Si = y~^ (see also Remark [T7)) . fv is not 
added to 7^2- 


Vi{x) 

^3- a 

Ui 


(i) r4: S2+j^ 


Sz+JK 


Vrr : 


S 2 


S 3 


(ii) yr^ = 82^ 83 = y'~ 


(iii) r((T3) = 


V 4 {x) 

a 

Ui^ 


■V 2 (x) 


V 4 (x) 

(I 

-U2, Ui 


U2 


V3(x) 


/ 0 The reac¬ 


tion 82 —S3 belongs to 7^2- 
(iv) Rate function: q{x){v2{x) + V3{x))v/i{x) 
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Hence, the reduced reaction network consists of two reactions: 




q{x)vi(x)vi{x) 


>S2 


q[x) (v2 (a;)+-!;3 (x))v4, (x) 


>83. 


Remark 19. If U is noninteracting, each reaction involves at most one species in U in the 
reactant and product. Each node Ui of a cycle a in Qu is the source (respectively, target) 
of exactly one edge of a. Therefore, in step ([4i|), each species in U that is a node of a cycle 
appears exactly once as reactant and once as product. 


Remark 20. We have defined a reaction network as a multidigraph to simplify the presenta¬ 
tion, in particular in connection with Theorem 1161 and hence also allowed multiple reactions 
between the same complexes. The reduced reaction network might, however, be further sim¬ 
plified by collapsing multiple reactions into one reaction without changing the ODE system, 
that is, using Q as explained in Remark [2j The rate function of a collapsed reaction is the 
sum of the rate functions being joined. See Example 1211 

Example 21. We continue with Example[9j The two components of the graph Gu are strongly 
connected, thus U is linearly eliminable. The set A consists of two cycles: ci, with the edges 
with labels kixs and k^, and a 2 , with the edges with labels k^xs and k^. Both cycles give rise 
to the reaction S —Sp and the rate functions are respectively 


Ki(x) = qi{x)kik 3 Xs and ^ 2 ( 3 ^) = q 2 {x)kikQXs, 
for qi{x) j[kiXs + k 2 + k^) and q 2 {x) {k^xs + k^ -\-kQ). By Remark [20l we might 

further simplify the network to S — Sp. 

Example 22. Consider the following reaction network: 


5i + Ui S 2 + U 2 S 3 + U 2 Si + Ui Si + U 2 S 3 + U 3 . 


Eor the noninteracting set U = {Ui, U 2 , U 3 } and a W-linear kinetics, the multidigraph Gu 

V2ix} 


has a spanning tree rooted at C/ 3 . Thus C/ is linearly eliminable. Consider the unique cycle 
O' of Gu- The graph a \ e with e: C/i —>■ C /2 cannot be extended to a spanning tree rooted at 
s(e) = Ui- Therefore, r((T) = 0 and the reduced reaction network obtained by elimination of 
U has no reactions. Since the graph Gu does not have a spanning tree rooted at a node other 
than C/ 3 , we have that ip{x) = (0,0, T), for the conservation law tti -|- ti 2 -|- tt 3 = T. Since the 
kinetics is CC-linear and C /3 is not involved in any reactant, it follows that k{x, (p{x)) = (0,0,0). 
Thus the ODE system of the reduced reaction network (jl6p is x = 0. 


5 Basic properties of the reduced reaction network 

In this section we study some basic properties of the reduced reaction network in Definition 
ESI In particular. Section [5.11 is concerned with the kinetics of the reduced reaction network. 
Section IQ is concerned with conservation laws. 

In what follows, a reaction network {C,TZ) on 5, a kinetics k and a linearly eliminable set 
U G S are given. We let (C, TZ) be the reduced reaction network on obtained by elimination 
ofZ^. 
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5.1 Rate functions and kinetics 


In this section we prove that if the original rate functions vanish whenever one of the concen¬ 
trations of a reactant is zero, then so do the rate functions of the reduced reaction network 
(Proposition [23]). A refined result is obtained when the kinetics of the original reaction net¬ 
work is mass-action (Proposition [25|). The proofs of Propositions 12311241 and 1251 are given in 
Subsection 19.41 

As noted in Remark [3l a biochemically meaningful rate function fulfils Kr{x) = 0 whenever 
the concentration of one of the reactant species is zero, that is, supp(yr) 2 supp(x) implies 
Kr{x) = 0. Based on this, we say that a function / defined on a subset V C M"" is standard 
for r G 7^ if 

for all X G V: supp(yr) ^ supp(x) f{x) = 0. (19) 

We say that the function is fully standard if the reverse implication in (11911 holds as well . 
The rate function Kr is (fully) standard if it is a (fully) standard function for r. Similarly, the 
kinetics k is (fully) standard if the rate function of all r G 7^ is (fully) standard. 

The next proposition states that if k{x,u) is standard, then so is k{x). 

Proposition 23 (Standard kinetics). (i) Let a G A and consider the reaetion r„ G 7 ^2- If 
the kinetics k{x,u) is standard (resp. fully standard) on Q x then n^a) is standard 
(resp. fully standard) on for r^- 

(a) If the kinetics k{x,u) is standard, then so is k{x). 

The kinetics of the reduced reaction network might not be fully standard, even if the 
kinetics of the original network is fully standard. For a reaction r G TZi defined by r G lZ\Itu) 
the rate function Kf(x) is fully standard provided Kr{x,u) is. However, other things might 
occur for reactions rv G IZ 2 , where the rate function might vanish if 

7 r(7 \ cj) = 0 for some x G Ll, 

7 er((T) 

cf. equation (|18l) . We characterise when this happens in Proposition 1241 below. Let rv be 
a reaction and x G M^o be positive in the coordinates corresponding to the reactants in 
TV. We show below that the rate function of To- vanishes at x, if and only if, the steady state 
concentrations of the species in the cycle a are zero whenever the concentrations of the species 
in are given by x. 

Proposition 24 (Fully standard kinetics). Assume that k{x, u) is fully standard onLlx M>q. 
Let a G A and x G Ll be such that supp(?/f^) C supp(x). Then the following two statements 
are equivalent: 

(i) n(iT) vanishes at x. 

(a) y^j{x) = 0, for all j G {1,, m} such that Uj gU is a node in a. 

It follows from the proposition that the rate function of To- vanishes if and only if ^^(e) {x, (/^(x)) 
0 for all e G CT with s(e) = Uj G lA. Indeed, the kinetics is 77-linear and hence Uj is a factor of 
the rate function of the reaction r(e). 

If the vectors x G LI for which Proposition [2H applies are precisely those that satisfy 
Xj = 0 for a certain i G {1,... ,p}, then the reaction To- can be modihed by adding Si to its 
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reactant and product, such that becomes fully standard for the modified reaction. This 
modification does not alter the ODE system. We do not apply this modification here. 

If the original reaction network is endowed with mass-action kinetics, then we have a 
sharper result than Proposition [23l In particular, the function vr(cr) has the form of a mass- 
action rate function for rv. Whether or not the rate function of the reaction is of 

mass-action type depends on the specific form of and ^ 7r(7 \ a). 

7 er((T) 

Proposition 25 (Mass-action kinetics). Assume that the reaction network {C,TZ) is endowed 
with mass-action kinetics. For a G A, 

7r{a) = X gQ, 


for some constant k„ > 0. Moreover, if the reactant of r^ involves the species Si, then the 
concentration Xi is a factor ofIl{cr) with exponent at least {yr,f)i- 

The proposition says that the exponents of the concentrations in n((T) might be larger than 
their stoichiometric coefficients, for the species involved in the reactant of r^. There might 
also be factors that depend on concentrations of species that are not involved in the reactant, 
as illustrated in the next example. The example illustrates also the following phenomenon: 
when two identical reactions arising from different cycles in Qu are joined and their respective 
rate functions are added (cf. Remark 1201 ) . then the rate function of the reaction might become 
fully standard. 

Example 26. Consider the reaction network with mass-action kinetics: 

Si + Ui^S2 + U2 U2 S3 + U3%Si + U2 U2 ^ U3 ( 20 ) 


and the noninteracting set U = {Ui, U 2 , U 3 }. The multidigraph Qu is 


k\x\ 



k2 k^xs 


The reduced reaction network obtained by elimination of Id is 

5i 52 53 5i (21) 

where k(x) = (ki(x), K 2 (®)) is 

k{x) = q{x)kik3(k2XiX3,kiXiX3), q{x) = —- ^ -- 7-7 -■ 

V / k2k3X3kik3XiX3kik4,xi 

Concentrations that appear with an exponent larger than the stoichiometric coefficient of the 
corresponding species are marked in bold. 

We have Q = i^>o \ {x I xi = X 3 = 0}. The rate function of the first reaction in ()21l) 
vanishes for X 3 = 0, even though S 3 is not in the reactant. This reaction corresponds to the 
left-cycle of Qu, which involves Ui,U 2 . By Proposition [2^1 the functions (pi{x),ip 2 {x) must 
vanish when X 3 = 0. Indeed: 

ip{x) =q{x) {k 2 k 3 X 3 , kiksxixs). 

Consider now the reaction network (|20p with an additional reaction: 

5i + 1/3 ^ 52 + U 3 . 


17 





The set U = {?7i, t/ 2 , C^s} is also noninteracting. The new reaction defines a self-edge at 
in Qu with label ^ 5 X 1 , and hence the reactions in (| 21 ll are also reactions of the new reduced 
reaction network with the same kinetics. The cycle given by the self-edge belongs to A and 
gives rise to an additional reaction —> S 2 with rate function q{x)kikjif,x\. This reaction is 

considered now twice with different rate functions. After joining the two reactions Si —)■ S 2 , 
the reduced reaction network is m with the rate function ki(x) replaced by: 

Ki(x) = q{x)kik2k3XiXs + q{x)kik^k^x\ = q{x)ki{k2k^x^ -|- k4k^xi)xi. 

With this transformation, the kinetics of (|21ll is fully standard in n. 

5.2 Conservation laws of the two networks 

Let S' C RP be the stoichiometric subspace of the reduced reaction network. In this section we 
compare the orthogonal complement of the stoichiometric subspaces, S-*- and S-*-, which define 
the conservation laws of the reaction networks. We show in Theorem EZl that the projection 
of S-*- on is contained in S-*-, and that the two spaces agree if the connected components 
of Qu are strongly connected. For example if xi -|- X 2 -k X 3 = T is a conservation law of the 
original reaction network and S 3 is eliminated, then xi -|- X 2 = T is a conservation law of the 
reduced reaction network. 

Let W be the vector subspace of S-*- generated by the vectors defined in (fTOjl for each 
connected component Qu of Qu such that * ^ M-u- Any vector in S'*- with support in lA belongs 
to W. That is, W is the kernel of the map Q : S'*- —RP (the restriction of the projection in 
dHI) to S'*-). The proof of the next theorem is given in Subsection 19.41 

Theorem 27 (Conservation laws). (i) C{S^) C S-*-. 

(a) If every connected component of Qu is strongly connected, then C{S'^) = S'*-. 

Theorem EH gives that C induces an injective map from the quotient vector space S-^/W 
to S'*", which is an isomorphism when all connected components of Qu are strongly connected. 
In this case, a minimal set of conservation laws of the reduced reaction network is simply 
obtained from a minimal set of conservation laws of the original reaction network by “deleting” 
Ui, i = 1 ,... ,m, from each conservation law, and any conservation law only in u. 

Example 28. By deleting ui = X 4 and U 2 = X 5 from the first equation in the minimal set 
of conservation laws ([7]) in Example 01 we obtain the minimal set of conservation laws of the 
reduced reaction network in Example 1181 that consists of 

Xl+X2 + X^= T. 

Example 29. We consider a mechanism that consists of two substrates A, B that are con¬ 
verted into two products P, Q through a series of reactions catalysed by an enzyme E [21 
Chapter 5]. It is an example of a bi-bi model in the notation of Cleland [ 6 ]: 

E + A^EA EA + B ^EAB ^EPQ ^EQ + P EQ^E + Q. 

k2 k4 ke ks kio 

We let f7i = E, U 2 = EA, Us = EAB, U 4 = EPQ and E^ = EQ. We assume mass-action 
kinetics. A minimal set of conservation laws is 

Te = Ui + U2 + Us + U4 + U5 Ta+q = xa + xq + U2 + Us + U4 + Us 

Tb+q = xb + xq + Us + U 4 + Us Tb+p = xb + xp + U 4 + Us. 


18 




We consider the reduced reaction network obtained by elimination of the set lA = {[/i, C/ 2 j 1 / 4 , C/ 5 }. 
The multidigraph Qu is strongly connected 


kiXA ksXB fcs kr 

Ui U 2 t /3 [/4 t/ 5 . 



By Theorem 1271 a minimal set of conservation laws for the reduced reaction network is 
xa + xq = Ta+q, xb + xq = Tb+q-, xb + xp = Tb+p- 

Let us hnd the reduced reaction network. A contains two cycles: the cycle with the edges 
with labels kiXA, k^XB, k^, kj and /cg that meets all nodes clockwise, and the cycle with 
the edges with the rest of the labels that meets all nodes anti-clockwise. For both cycles, 
r(iT) = {fj} (Remark [2]) and the reduced reaction network is: 


A + B 


K2(a:) 


P + Q 


where ki{x) = q{x)kik 2 ,k^kjk^XAXB-, i^ 2 {x) = q{x)k 2 k 4 kQkskioXpXQ and q{x) equals Tp di¬ 
vided by 


{k/^k^ + kik7 + k^kj){k2kQ + kik^XA + k2kioXQ) + kikQksXp{kiXA + ^ 2 ) 

+ ksk^kYXpikg + kioxq) + {kjik^ feg) -b (/cs + kQ){kg -b k^xp))kik2,x axb 
+ {ki{k2 + fee) + (fes + fe6)(fe2 + k^XB))k%kiQXpXQ. 

In this case = M>g. The reduced reaction networks for this example and for the ping-pong 
bi-hi mechanism considered in the introduction have the same reactions. The two networks 
differ only in the factor q{x). 

If some connected components of Gu are not strongly connected, then the cases S'-*- = C('S'-’-) 
and S'-*- 7 ^ C('S'-*-) are both possible. 

Example 30. Consider the reaction network in Example 1221 A minimal set of conservation 
laws consists of 

U 1 +U 2 + U 3 = Ti, Xi-\-X 2 = T 2 , + 3^4 = Ls, Xi - X 3 - Ui U 3 = T 4 . 

By Theorem 1271 the following are conservation laws of the reduced reaction network: 

3^1+3^2 = 12 , 3:3-bX4 = r 3 , Xi - X 3 = r 4 . 

However, the reduced reaction network has no reactions and therefore S'-*- = M^. Here, the 
multidigraph Gu is not strongly connected and C{P'^) $ <5-*-. 
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Example 31. Consider the reaction network 

s, + S,_ + c/2 S2 Si 

and the set U = {Ui,U 2 }- A minimal set of conservation laws of the reaction network consists 
of ui + tt 2 = Ti and xi + X 2 = T 2 . The multidigraph Qk is 

c/i "W > C/2 

which is acyclic. The reduced reaction network has one reaction S 2 —^ /S’!, and a minimal set 
of conservation laws consists of xi + X 2 = T 2 . The multidigraph Qu is not strongly connected 
and S'-*- = C(‘S'''‘)- 

6 Iterative elimination 

In this section we discuss stepwise elimination and show that in certain cases one might obtain 
different reduced reaction networks, depending on the order by which species are eliminated. 
This naturally has practical implications. Biologically, it relates to the situation in which 
species at different time-scales are eliminated from the network. The fastest species are re¬ 
moved first, then the second fastest and so on. As the resulting network might differ from 
the network obtained by removing the species all at once, a proper partitioning of species into 
time-scales might be warranted. 

The setting is the following: Let 

0 = Wo C ZYi C ZYs C • • • C Wz = W 

be an increasing sequence of linearly eliminable sets. We compare the situation where W is 
eliminated at once (direct elimination) to the situation where the sets Wj \ Wj_i, i = 1 ,... , 1 , 
are eliminated iteratively (iterative elimination). 

If each set Wj \ Wj_i, i = is the node set of a connected component of Gu (or a 

union of such sets), then direct and iterative elimination yield the same result by construction: 
both elimination of variables and the reduced reaction network consider the multidigraph Gu 
componentwise. 

Another scenario in which iterative and direct elimination coincide is as follows. Consider 
the connected component Gh of Gu that contains * (if any) and let be the multidigraph 
obtained by removal of the node * and all the edges that have * as source or target. This 
multidigraph might now have several connected components. The next proposition (proven 
in Subsection 19.5p states that elimination of the species in "H or iterative elimination of the 
species in the connected components of G^ yield the same reduced reaction network. 

Proposition 32 (Splitting the connected component with node *). Assume that Gu is con¬ 
nected and contains *. Let U = LLi U 7^2 a decomposition of U into two disjoint subsets 
such that Gh = Gfii U Gfi^ ■ Then the reduced reaction network obtained by elimination of U 
agrees with the reduced reaction network obtained by elimination ofLLi followed by elimination 

of 712. 

See Example [3U] in Subsection 17.21 for an illustration. 

In general, iterative elimination might not yield the same reaction network as direct elim¬ 
ination. We present some examples below. 
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Example 33. Consider the reaction network with mass-action kinetics 


Si^Ui + S 2 ^U 2 ^Ui+S 3 % S4, 


and the linearly eliminable sets Ui = {Ui} CU 2 = {Ui, U 2 } with multidigraphs 

Gui ■ Gu 2 ■ 



koxo 

Ui^U 2 


\i 


The four cycles of Gu^ give that the reduced reaction network obtained by elimination of Ui is 

( 22 ) 


S, + S2 S2 + C/2 S3 + U2 Ss + St 


S, + s,^i^S2 + St 


with kinetics 


^ 2 _ f kik2XiX2 k2k^X2U2 ^ 3 / 04 X 3^2 kiktiXiX'i \ 

\k 2 X 2 + kiX2, ’ k2X2 + kt^X^ ’ k2X2 + kiXz ’ k2X2 + kiXs J ' 

Here H = \ {x 2 = X 3 = 0}. The set W 2 \^i = {^^ 2 } is linearly eliminable for the reaction 


network ([2^ . Its elimination yields the reduced reaction network: 

Kn(x) 

S 2 ^3 + ^3 


^1+52 + ^3 ^2 + ^3 + ^4 


K'g{x) 


>52 + 54 


with rate functions 


k\{x) = 

k^ix) = 


kik2k^kiXiX2X^/{k2X2 + ^ 4 X 3 )^ kik2XiX2 


k^k/^xi,/{k2X2 + ^ 43 ^ 3 ) k2X2 + k^xs ’ 

k2k3X2/{k2X2 + k4^X3) _ ^ 2 X 2 
^ 3 ^ 43 ^ 3 / {k2X2 + A:4X3) kiXs 


~/ / \ ^ 4 / 04 x 4 x 3 ~ ~ 

In this case H' = {x G M>q | X 3 7 ^ 0}. Since ^ 2 X 2 + k^x^ 7 ^ 0 for all x € H' we might simplify 

k 2 X 2 + A: 4 X 3 from K 4 (x) and ^^(x), as shown above. 

The reduced reaction network obtained by direct elimination of U 2 is: 


with rate functions 


Kl(X = 


S 2 ^S 3 

54 + 53 ^ 52 + 54 

^2^3X2 /C2X2 

^3/04X3 ^4X3 ’ 

_ ^4/04x4x3 

K2(x) = 

/C3/C4X3 


k3 ’ 


which is clearly different from the reduced reaction network obtained by first eliminating Ui 
and then U 2 \Ui. 
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When Qu is connected and contains elimination of u only involves the equation u = 0 
and no conservation laws. Therefore, iterative elimination is equivalent to solving a system 
of linear equations iteratively. It follows that the function y? in (fT^ does not depend on the 
chosen procedure and thus neither does the ODE system ()16h . However, as shown in the 
above example, the reduced reaction network may depend on the procedure, even though 
their associated ODE systems agree. 

The next example provides an example where direct and iterative elimination result in the 
same reduced reaction network, and an example where the reactions but not the kinetics of 
the two reduced reaction networks agree. 

Example 34. Consider the reaction network with mass-action kinetics: 

Si + Ui%U 2 %U 3 % S2 + Ui, 

and the sets of noninteracting species: 

Ui = {Ui}ciU2 = {Hi, U2} C 1 U^ = {C/i, C/2, C/3}. 

The multidigraphs Qui are 



The reduced reaction network obtained by elimination of lAi is 

5, + [J3 S2 + f/2 t/j U, ( 23 ) 


with rate functions 


Ki{x,U2,U3) 


hksxius 

kixi 


ksus, 


K2{x,U2,U3) = k2U2. 


Since D = {x € M>q | xi 7 / 0}, xi can be canceled in ki(x, U 2 , ■us). 

The set C /2 \ C^i = {C/ 2 } is linearly eliminable for the reaction network (1231) . The multidi¬ 
graph Qu 2 \u^ is 





Elimination of IA 2 \ IA\ yields the reduced reaction network 


5i + C/3 


K.[{X,U3) 


> ^2 + C/3, 


k[{x,U 3 ) 


k2k3U3 

k 2 


k3U3- 


(24) 


This is also the reduced reaction network obtained by elimination of C /2 directly, as it can 
easily be seen by considering Therefore, in this case both eliminations yield the same 
result. 

We compare now direct elimination of C/3 with elimination of C/3 \ C/2 = {C/3} from (|24p . 
Both approaches provide the reduced reaction network 




k ' i ' (x) 


> 52 . 
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However, the kinetics differ: k"{x) = k^T, for iterative elimination, using the conservation law 
U 3 = T, and 


kik2k‘iTxi 

kik2Xi + k2X3 + kik^xi ’ 


using U 1 + U 2 + U 3 = T, 


for direct elimination. Therefore, the kinetics of the reduced reaction network depends on 
whether the elimination is performed iteratively or not. 


In general, when conservation laws are involved in the elimination procedure, then the 
kinetics obtained after direct and iterative elimination differ, even though the reactions of 
the two reduced reaction networks might be the same. The reason is that the reduced ODE 
systems x = C{x,Lp{x)) in (1161) obtained by iterative or direct elimination are different. 

To understand why, let U = {Ui,... ,Um} be a linearly eliminable set such that Qu is 
strongly connected and does not contain the node *. Then ui + • • • + Um is conserved. The 
set Ui = {Ui,... ,Um-i} is linearly eliminable and Gu^ contains the node *. Theorem [TT] 
thus guarantees that the system of equations iii = ■ ■ ■ = iimi = 0 in ui,... ,Um-i has a 
unique solution ip^{x,Um)- The reduced reaction network obtained by elimination of has, 
by Theorem 1271 the conservation law 

Um — T. 

This implies that in the ODE system x = (^{x,ip{x)) obtained by iterative elimination of the 
set Ui followed by {Um}-, ^m{x) is constant equal to T, and the first m — 1 components of (p 
agree with ip^{x,T). 

If U is directly eliminated, then ip{x) is found by solving the system of equations 


ill — • • • — iim—i — 0 , 


i=l 


The first m — 1 equations give Ui = ipj{x,Um) for z = 1,... ,m — 1, and the last equation 
determines Um as the solution to 


m —1 

P\{X, Um) + Um=T. 

i=l 


m —1 

Unless is constant, Um = P>m{x) is not constant. Hence the function ip{x) used 

i=l 

in i; = (^{x^ip{x)) depends in general on the chosen approach, and as a consequence so does 
the kinetics of the reduced reaction network. 

7 Relation to previous work 

Ideas similar ideas to ours have been proposed in the literature. In the first part of this section 
we focus on the techniques proposed in UnillZ]; see also [mug]. We will discuss similarities 
and differences to our approach following the exposition in [18]. In the second part of the 
section, we justify that our construction generalises the reduction of intermediates in 1 ^. 
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7.1 Horiuti-Temkin approach 

In [IB] the authors outline some applications of graph theory to the theory of reaction networks 
and introduce reduction of reaction networks by graphical means. In that work, the reversible 
reactions are treated as one reaction, while we treat them as two separate reactions in the 
present paper. 

The authors introduce a graph called the kinetic graph (see m) whose nodes are the 
so called intermediates. The kinetic graph coincides with the multidigraph Qu introduced in 
Definition [7] (up to the treatment of reversible reactions) and the so called intermediates form a 
noninteracting set in our terminology. Their goal is to eliminate intermediates from a reaction 
network and find a minimal mechanism that allows the computation of the production rates 
of the remaining species. The differences between their approach and our work arise from the 
details in the treatment of the reactions and the rates. 

First, in [TB] not only one reduced reaction network is obtained, as in our case, but an 
infinite number of them. Any linear combination of the original reactions that cancel the 
intermediates is a possible reaction in a reduced reaction network. The set of such reactions 
defines a vector subspace of S, of which a basis is chosen. Therefore, the reduced reactions 
are independent vectors in the stoichiometric subspace, and their number is minimal. 

Second, the conditions imposed for the computation of the rate functions are different. In 
|18j . the rate functions are found by imposing the Horiuti-Temkin equation (see (I25p i. which 
involves the rate functions of the original and reduced reaction networks but does not relate 
to the stoichiometry of the reactions. As a consequence of the Houriti-Temkin equation, the 
ODE system (fTB|l is also satished for the reactions and rate functions given in [TB] . We discuss 
next in further detail the differences on this particular point. Assume for the discussion below 
that TZ = TZy. 

For fj € A, let be the vector with = 1 if Cj is an edge of the cycle and (t'o-)i = 
0 otherwise. Choose an order for the set of edges in Qu and for the set of cycles A = 
{fji,... ,cr|/^|}, and let H be the |A| x i matrix whose ith row is With this notation, 
the Horiuti-Temkin equation for the reduced reaction network in Definition [T5l considering 
reversible reactions as two irreversible reactions, reads 

k{x,(p{x)) = H^k{x). (25) 

Componentwise, this condition is in our notation 

Ki{x,(p{x))= = X] z = 

aeA,e(ri)£(T 

This condition might not be satisfied by the rate functions in Definition [T5j We show this in 
the case where Id is linearly eliminable and Gu is connected and does not contain *. We let 
q{x) = qu{x). 

Let z G { 1 ,... ,£} and Uj € Z7 be involved in the reactant of rj G 77 (it exists because Gu 
is connected and does not contain *). Then 

Ki{x,(p{x)) = ipj{x)vi{x) = q{x) ^ 7r(r)z;i(x) = (*). 

Tee(Uj) 


Each term 7r(r)z;j(x) in the sum is the label of an element in F(f7) for some cycle a that 
contains the edge e(rj). Let A' be the set of all cycles of Gu- Using the definition of F(it) we 
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have 


w = q{x) Y1 Y1 

uGA' ,r{iT)f^@,e(ri)€cr crG A,e(ri)Gcr 

The inequality arises because Y1 Civ'rie) ~ Vrie)) 7 ^ 0 for cr G A and thus the sum on the 

left-hand side might involve more terms than the sum on the right-hand side. Hence, the 
Horiuti-Temkin equation is not necessarily satisfied. 

Our rate functions fulfil an equation similar to the Horiuti-Temkin equation, once stoi¬ 
chiometry is introduced. To understand this, write the ODE system ([ 6 ]) as 


g{x) = Ak{x), X G H, (26) 

where A = {aij) G is the stoichiometric matrix with = {y'^. — yrj)i- The stoichiometric 
matrix of the reduced reaction network is by Definition [15] where A^ is the matrix given 

by the first p rows of A. By equation (12611 and Theorem fTHl the kinetics in Definition [15] 
satisfies 

(/?(x))) = A'^H^k{x). 

Finally, the kinetics obtained in [18] is not necessarily standard if the kinetics of the original 
reaction network is, contrary to our kinetics (cf. Section l5.ll) . This is illustrated using the 
main example in [18] . 

Example 35. Consider the reaction network with mass-action kinetics 

Si + S2 + Se^ S3 + S7 S7 ^ Si + Se 84 + 87 ^ Si+ S5 + Sq 

k2 ^4 ^6 

The correspondence with [TH] is as follows: 5i = 0,^2 = H 2 O, S 3 = H 2 , ^4 = CO, S 3 = CO 2 , 
Sq = Zi and S7 = COZi. We consider the linearly eliminable set U = { 56 , 5 ' 7 } with the 
conservation law xq + X 7 = T. The reduced reaction network obtained by our procedure is 


Kl{x) 

Si + S2 S3 + Si 

K2{x) 


Si + 82 + Si Si + S3 + S3 

K 4 ,{x) 


2 Si 


K+X)^ 
\ - 

K+X) 


S 1 +S 3 


with kinetics 

k{x) =q{x) {kik 3 XiX 2 ., k 2 kiX 3 Xi, kik 3 XiX 2 Xi, ^ 2 / 06 X 1 X 3 X 5 , kik 3 xl, /cs/ce^ixs^, 


where q{x) = T{kiXiX 2 + kQXiX 3 -|- k 2 X 3 + {ki -|- A; 5 )x 4 -|- /cs) ^. The kinetics is standard for the 
reduced reaction network, because each rate function vanishes when one of the concentrations 
of the species in the reactant is zero. 

One of the reduced reaction networks obtained in [T 8 | is 


S1 + S2 


k[{x)^ 

T- 


^3 + ^4 


52 + ^4 


k'+x)^ 

T- 

k'^{x) 


S3 + S3. 


(27) 


In [T 8 | the rates of the reactions are considered by pairs of reversible reactions and we obtain 
by their algorithm that 

“ ^2(^) =Q{x){kik3XiX2 — k2kiX3Xi -|- /cs/oexixs — kik3Xi) 

K 3 (x) - K'i{x) =q{x){kik 3 XiX 2 X 4 - k2k3XiX3X3 - k3k3XiX3 -I- kik3xl) 


with q{x) as above. By collecting the terms according to their signs, we find 

Ik'iix) =q{x){kik 3 XiX 2 + k 3 k 3 XiX 3 ), 'k' 2 {x) =q{x){k 2 kiX 3 X 4 + kik 3 xl). 


This kinetics is not standard for the reaction network (|27p . 
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7.2 Intermediates 


Our construction generalises the reduction of intermediates in [3]. In [3], an intermediate 
is defined as a species y in a reaction network that is produced in at least one reaction, 
consumed in at least one reaction and is not involved in any complex other than Y. A set 
of intermediates y is a subset of the species set and at the same time a subset of the set of 
complexes (under the identification of complexes with linear combinations of species). 

Any set of intermediates is noninteracting. The multidigraph Qy contains the node * and 
is strongly connected. Thus, assuming mass-action kinetics, T is a linearly eliminable set. 
The reduced reaction network obtained by elimination of y has the following properties. 

Any cycle u € A must contain the node because at least one reaction r(e) with e in cr 
must involve a species in y^. Then the reactant (resp. product) of To- is the reactant y^^e) 
(resp. product of the reaction corresponding to the edge e with source s(e) = * (resp. 

target t{e) = *). 

The label of an edge in Qy whose source is not * is constant. Thus if r is a spanning tree 
rooted at *, then 7r(r) is a product of reaction rate constants /cj. If r is a spanning tree rooted 
at y € y, then 7r(r) is a product of reaction rate constants and the label of the edge e 

with s(e) = * in the tree. Thus 

q{x) = G{k)~^, and II((t) = 

where G{k),F{k) are polynomials in the reaction rate constants, such that q{x) is constant in 
X. Therefore, the reduced reaction network obtained by elimination of y is a reaction network 
with mass-action kinetics. There is a reaction between two complexes in C \ y in the reduced 
reaction network if and only if the reaction is already in the original network or there is a 
directed path between the two complexes through intermediates. 

Example 36. Consider the reaction network with mass-action kinetics: 

+ 52 ^ Tl ^ Ta + Sg Tg ^ Tg ^ ^4 + 56 

k2 ks 

+ S4 + 55 5 g + 57 ^ 52 + 57. 


A set of intermediates of this network is y = {Ti,... ,y 5 }. The corresponding multidigraph 

Gy is 



By Proposition 1321 the reduced reaction network obtained by elimination of y can be found 
by iteratively eliminating the sets of nodes yi = {yi,y 2 } and y 2 = {yg,l 4 ,y 5 }. The multidi¬ 
graphs Qy^, Qy^ correspond to the left-subgraph and right-subgraph of Qy^ respectively. We 
obtain the following reduced reaction network with the specified rate functions: 

5 1 + 52 ^ 5 i + 5 g 52 + 54 ^ 54 + 55 

52 + 54 ^ 54 + 56 5 g + 57 ^ 52 + 57 
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where 


Ki{x) 

K3{x) 


kik3kiXiX2 kik3XiX2 
ki{k2 + fca) (^2 + kz) ’ 
kekrksx^xr k^k-jx^x^ 

fc8fcio(fc7 + fc9) kio{k-j + kQ) 


_ kekgkipx^xr _ kgkgx^xr 

^ kskio{kr + kg) ksikr + kg)' 
K4{x) = kiiXgXT- 


There is indeed a reaction between every pair of complexes in C \ T that are connected by a 
directed path through intermediates. Further, the reduced reaction network has mass-action 
kinetics as well. 


8 Post-translational modification networks (PTMs) 

We conclude by discussing reduction of PTMs. In this section we abuse notation and use x 
for the vector of concentrations of the original species set as well as for the concentrations in 
the reduced reaction network. 

A common feature of signalling systems is the incorporation of PTMs, the attachment of 
some chemical group to a protein, after it has been translated. The most common example 
is phosphorylation. A mathematical formalism to study PTM networks, that is, a network 
combining several PTMs, was introduced in |20] . 

In [Soils] PTMs are considered from the point of view of variable elimination, where so- 
called substrates and intermediates are eliminated. This provides a system of equations that 
depends on the enzyme concentrations only. Here we study the reduced reaction network on 
the set of substrates obtained by elimination of the sets of enzymes and intermediates. 

We start by giving the definition of a PTM network, which is slightly more general than 
the one in [20115]. The species set of a PTM network is the disjoint union of three non-empty 
species sets: 

• a set of substrates S = {S'!,..., S'p}, 

• a set of enzymes £ = {Ei ,..., T'mi} and 

• a set of intermediates y = {Yi ,..., Tmj}, in the sense of Subsection 17.21 
Allowed reactions, taken with mass-action kinetics, are of these five types: 

+ Ej —^ Yi, Yi—^S^ + Ej, Y, —^ Yj, 

Si + Ej Si + Ej. 

The bottom types are not considered in We assume that any path 

5,, + Ej, - S,, + Ej, (28) 

through intermediates satisfies ji = j 2 - This provides a decomposition of y into at most mi 
disjoint subsets according to the enzyme that ultimately produces them, or to which they 
dissociate. 

With lA = £\jy, the multidigraph Qn has a connected component for each enzyme, which 
is strongly connected by the hypothesis and does not contain *. Thus lA is linearly eliminable. 

The next proposition (proven in Subsection 19.6p describes the reactions of the reduced 
reaction network {C,TZ) obtained by elimination of lA. Recall the decomposition TZ = TZi U7^2 
given in Definition [TSl A reaction Si —Sj in the reduced reaction network belongs to TZi if 
it already belongs to TZ. Therefore, we only consider 7^2 in the proposition below. 
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Proposition 37 (PTM networks). Let {C,TZ) be a PTM network and {C,TZ) be the reduced 
reaction network obtained by elimination ofU = Euy. Let & S be two substrates. 

The reaction Si^ —>■ Si^ belongs to 7^2 if and only if there is a path as (1^51) from Si^ + E to 
Si 2 + E in (C, TZ) for some E ^ S. 

We conclude that if two substrates do not interact with a common enzyme E, there is no 
reaction between them in 7^2- 

Example 38. Consider the network in Example [36l which is a PTM network with y as in 
the example, £ = {^i, 5'4, S'y} a set of enzymes and S = {S 2 , S 3 , S^, Sq} a set of substrates. 
By Proposition 1371 the reactions of the reduced reaction network obtained by elimination of 
£uy are 

^2 ^ ^3 ^2 —^ 55 ^2 —> ^ 6 . 

Example 39 (n-site phosphorylation system). We consider an n-site sequential distributive 
phosphorylation mechanism, which consists of a substrate S that contains n ordered phospho¬ 
rylation sites. We let So denote the unphosphorylated form and Si denote the phosphorylated 
form in which sites 1 to i are phosphorylated. We assume there is a kinase E that catalyses 
all phosphorylation steps and, similarly, a phosphatase F that catalyses dephosphorylation 
steps. The reaction network associated with this system is 


So + E 

Sn + F 


T- 

K,o(x) 

kI(x) 


kUx) fdAx) 

Yi Si+E ^2 — 

kI{x) 

K.i(x) kI{x) 

Zl + F Z 2 

kI(x) 





Sn y E 




5o + F. 


Let U = {E, F,Yi, ..., Y^, Zi,. .., Z„} and assume the kinetics is 77-linear. By Proposition [371 
the reduced reaction network obtained by elimination of lA is: 


k\(x) 

Sn^^Si 

kA{x) 


"LiP) 


nlix) 


Sr,.. 


In order to find the kinetics, we consider the multidigraph Qk, which has two connected 
components: 




The rate functions Vi depend on x. The connected components give rise to conservation laws 
with total amount and T^, respectively. The elements of A are: 


<j}- 




Zi 
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for i = 1,... , n. For I = 1,2 and i = 1, ..., n the reaction with rate function k\{x) corresponds 
to the cycle u* and is as follows: 

n 

k[{x) = (U3j_i + V^3j) 

n n n ^ 

qi{x) = + X!/ ^3i-2 JJ + ^Sj)) ■ 

j=l i=l j=lj^i 


9 Proofs 

In this section we present the technical details and proofs of the results stated in the previous 
sections. This section is organised such that it follows the structure of the rest of the paper. 


9.1 Preliminaries 

We give the precise definition of the digraph Q introduced in Section [2] and some associated 
functions that will be used later on. 


Definition 40. Let Q = {Af,£) be a labeled multidigraph and let £' = {e ^ £ \ s{e) ^ t(e)} 
be the set of edges that are not self-edges. The labeled digraph Q = {M, £) associated with Q 
is the graph with 

M = {A^ E AA I 3e ^ £' with t{e) = N or s{e) = N} 

£ = {W N 2 I 3e ^ £' with s(e) = Ni and t{e) = N 2 }. 


Note that by definition J\f C J\f and that the inclusion might be strict if there is a node in 
Q that is only connected to itself. We define a surjective map from £' to £ as follows: 


£' ^ £ 

e !-)• s(e) —>■ t{e). 


(29) 


Using this map, the labeling tt for Q is defined as 

eG/3“i(U 


9.2 Elimination of variables 

Theorem m is stated and proven in [3] using the digraph Qn instead of the multidigraph Qu 
as we do here. This affects the definition of q'u{x) as well as the sets of spanning trees. Either 
way the functions (p agree because the computations performed using Qu or Qk agree. It is 
shown in the next lemma. 

Lemma 1. Let Q = [M,£) he a labeled multidigraph and let Q = {M,£) he the associated 
digraph given in Definition Let N £ Af and let Q{N), B{N) be the set of spanning trees 
rooted at N of Q and Q, respectively. Then 

^ ^{ t ) = 7r(r). 

Eee(v) T&e{N) 
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Proof. Since a spanning tree cannot contain a self-edge, the map {3 in (IMl) extends to a 
surjective map from Q{N) to Q{N). In particular, 


0(iv)= ]j rHn 

fG0(iV) 

Let r be a spanning tree rooted at N in Q. Then, 

^f(r) = n ^(e) = n E ^(e) = E 0 7r(e) = ^(^)- 

e£T eG/3“l(e) rE/3“^(r) 

Therefore, using (l30|) we obtain 

rG©(7V) Tee{N) 'rG/3-l(^) fe©(Af) 


(30) 


□ 


9.3 Reduced Network 

In this section we first prove that the definition of r((T) in (I17p . 

r((T) = {rUe I a\e C T, r€ ©(s(e))}, 

is independent of the chosen edge e, cf. Proposition 1411 Subsequently, we prove Theorem 1161 
We define the following set of sub-multidigraphs of Gu- 

r := {7 = r U e I r G 0(s(e)), e^Gu}- (31) 

Each element of T is the union of an edge e of Gu and a spanning tree (of the connected com¬ 
ponent in Gu that contains e) that is rooted at the source of the edge, s(e). Any multidigraph 
7 G r contains a unique cycle. Indeed, for 7 = rU e as in ([STI) . the cycle is obtained by joining 
e and the path in r from t{e) to s(e), which exists because r is rooted at s(e). In particular, 
the cycle contains the edge e. 

Proposition 41. Let a be a cycle of Gu- For any edge e G a, 

{rUe I a\e C T, rG 0 (s(e))} = {7 G T | cr C 7}. 

Proof. Without loss of generality we assume that Gu is connected. If 7 G T and u C 7, then 
7 = r U e with r G 0(s(e)) and e G a. Clearly r is a spanning tree that contains a \ e. 
Conversely, if r G 0(s(e)) is a spanning tree containing a\e, then rUe G T and a C rUe. □ 

Therefore, r((T) is the subset of T whose elements contain a. This shows that the definition 
of r(cr) in (jl7l] is independent of the choice of edge e. Note that T is the disjoint union of 
r(cr) for all cycles a. 

The next proof shows that the reduced ODE system (fTHl) is the ODE system associated 
with the reduced reaction network in Definition [T5j 
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Proof of Theorem f7^ - Reduced reaction network. Let 

f{x) = ^r{x){y~ -yr)+ ^Ax){y~ - y^). 

rST^i re7^2 

We want to prove that f{x) = g{x) = C{g{x,(p{x))). Observe that 

9{x)= '^r{x,(p{x))C{y'r -yr) + ^ Kr{x, ip{x))C{y'r - yr) ■ 

r^TZu reTZu 

By definition of TZi and K^ix) for r € TZi, the first summand of f{x) and that of g{x) agree. 
Using the definition of 7 ^ 2 ; ah we need to prove is that 

'^QH{a){x)n{a)'^C{y'^^^^-yr(e)) = Hrix,ip{x))Ciy'r-yr)- (32) 

o-eA eGo- r&TZu 

The sums on the right- and left-hand sides can be decomposed into sums over the connected 
components of Qu that contain a or e(r), respectively. Therefore, it is enough to show that 
([32]) holds when Qu is connected. We let q{x) = q'^i^^){x) = qu{x). 

Given r € TZu-, either s(e(r)) = Ui for some i if p{yr) 7 ^ 0 , or s(e(r)) = * if p(yr) = 0 . In 
the former case Ui is the only species in U involved in the reaction r. Since the kinetics k is 
7Y-linear, we have for r G TZu 


Kr{x, ip{x)) 


ipi{x)vr{x) if s(e(r)) = Ui 
Vr{x) if s(e(r)) = *. 


By the definition of q{x) in (fT^ and of ipi{x) in (fTbl) . we have 


Ur{x,(p{x)) = q{x) E 7r(r) I Vr{x), for all r € TZy. 

\rG0(s(e(r))) / 

Comparing this equality with (1321) . the statement follows if the following holds: 

X] Z] T^i.x)Vr{x)C,{y'^ - yr) = '^'^i(^)Ciy'r{e)-yr(e))- (33) 

rGlZu TGS{s{e{r))) ctGA eEcr 

We show the equality from right to left. Let A' be the set of cycles in the multidigraph Qu. 
Since either ^ Ciy'rie) ~ hrie)) = 0 or n(iT) = 0 for u € A' \ A, we obtain that 

Z] Z]n(f^)C(2/(.(e) -2/r(e)) = Z '^'^(■^)Ciyr{e) - Vrie)) 

ctGA eGcr ctGA' eGcr 

= Z Z Z ’^(7\e)vr(e)C(2/(.(e) - yr{e)) = (*)• 

CTGA^ eE(T 'y£r((T) 

Each digraph 7 \ e in the sum is a spanning tree rooted at s(e). There is a bijection between 
the set of triplets (a, e, 7) such that <7 G A', e G a, 7 € r((T) and the set of pairs (e, r) such that 
e is an edge of Qu and r is a spanning tree rooted at s(e). Using further that 7r(e) = Vr(e)(x) 
and the correspondence Q, we obtain: 

(*) = Z Z ’^('^)^(e)C(2/r(e) -yr(e)) 

eeSu TG0(s(e)) 

= Z Z ^MVr(x)C(yr - yr)- 

tGTZu TG©(s(e(r))) 

This shows that ()33p holds, which concludes the proof. □ 
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9.4 Basic properties of the reduced network 

In this section we prove the results about the kinetics and the conservation laws of the reduced 
reaction network in relation to the original reaction network. 


Kinetics. 

Proof of Proposition [151 - Standard kinetics, (i) By definition 

yr. ='^C{yr{e))- 

eda eda 

Hence, 'k{(j) = 0 if and only if Vr(e){x) = 0 for some e G a. Further, we have that supp(?/f^) = 
U supp(C(yr(e)))- Assume that the kinetics k is standard on H x K>q. For x G H, we have 

eGa ~ 


supp(?/f^) 2 supp(x) 3 e G cr: SUpp(C(yr(e))) 2 supp(x) 
=> 3 e G (t: u^(e)(x) = 0 = 0. 


This shows that Tr{a) is standard on H for Pa- If the kinetics is fully standard, then the reverse 
of the second implication holds, showing that '7r{a) is also fully standard for Fo-. 

(ii) For Fo- G 7^2 defined by ct G A we have Kr^{x) = (jr(x)n((T). By equation (fTHl) . is 
a factor of n((T) and hence is standard on H. 

For F = dVr) Civ'r) € with r G TZ\R,u-, we have Kr{x) = Kr{x, ip{x)) by Definition fT5l 
Since r ^ TZu, then supp{yr) = supp(C(?/r)) = supp(7/f). Since Kr is standard on D x M>q, 
then Kr{x,u) vanishes if x G D fulfils supp(yf) ^ supp(x). Since the denominators of qn{x) 
do not vanish for x G D, we have Ff(x) = 0 whenever supp(yjr) g supp(x) and x G D, showing 
that Ff(x) is standard on D. □ 

Proof of Proposition 1^41 Fully standard kinetics. Throughout the proof all functions are 
assumed to be evaluated at x. Assume that statement (ii) holds. Let e G u be an edge such 
that s(e) = Uj G U. Then 


n(cj)=7r(e) ^ vr(7 \ e) < 7r(e) ^ 7r(r) 

7er((T) TdOiUj) 


= 7r(e) 


Pjjx) 

Quipid' 


Since q:>j{x) = 0, n((T) = 0 and statement (i) holds. 

Assume now that statement (i) holds. Let Uj G U he a node in a, r a spanning tree 
rooted at Uj and e G a the edge with source Uj. We construct a new tree r as follows: for 
every Ui ^ Uj in u, replace the only edge with source Ui in the tree r by the edge in cr with 
source node Ui. The obtained subgraph r is also a spanning tree rooted at Uj that satisfies 
r U e G F((t). 

By assumption (i) and the definition of n(cT), 7r(r U e) = 0. Since k{x, u) is fully standard, 
7r(e') ^ 0 for all e' G a by Proposition [23l Therefore there must be an edge e of r that does 
not belong to a and such that 7r(e) = 0. By construction, this edge is also an edge of r and 
hence tt{t) = 0. 

This proves that it{t) = 0 for all r G Q{Uj), which implies that ipj{x) = 0 (cf. 

Hence statement (ii) holds. □ 
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Proof of Proposition \2^ Mass-action kinetics. The label of an edge of a is 

where fcj.(e) > 0 is the reaction rate constant of reaction r(e). Using the definition of in 

Definition 1151 we obtain 


7r((T) = 

e£a 






Hence Tr{a) has the claimed form with = 0 ^r{e)- The last statement is a consequence of 

edcr 

(HID. □ 


Conservation laws. We prove here Theoreml27[ Before that, we introduce some graphical 
constructions and a technical lemma necessary for the proof. 

Let Q = {M,£) be a multidigraph and assume that the sets M = {Ni,..., Nn} and 
£ = {ei,..., Cf} are ordered. We define the following objects: 

1. The incidence matrix Cg of Q is the n x i real matrix such that 

(I if W = t{ej) 7 ^ s{ej) 

{Cg)ij = S “1 T ^i = / t{ej) 

I 0 otherwise. 


2. The cycle space of Q is the kernel of the incidence matrix. 


Let A' be the set of cycles of Q. For a € A', the vector t'o- with (r'o-)i = 1 if Ci is an edge of 
the cycle and = 0 otherwise, belongs to the cycle space of Q. Moreover, the elements t'o- 
correspond to the irreducible elements in the cycle space with all non-zero components equal 
to one. That is, cannot be expressed as the positive sum of two vectors in the cycle space 
with non-negative integer coordinates. The elements in A' are also called elementary cycles 
in the literature [22]. We choose an order for the set of cycles A' = {cji,..., cr|^/|}, and let H 
be the |A'| x i matrix whose ith. row is 


Lemma 2. Let Q be a strongly connected multidigraph. Then kerLf = imCg. 

Proof. By duality it is enough to show that imH* = kerCg. Since for all a € A', z/g- belongs 
to the cycle space, we have CgX = 0, that is, imX C kerCg. Since Q is strongly connected, 
then there exists a basis of ker Cg whose elements are of the form for u € A' (see [T] ). 
Hence im X D ker Cg as desired. □ 

We have now all necessary tools to prove the relation between and S-^. 


Proof of Theorem 21: Conservation laws. (jl|) Let oj € and r £ IZi be a reaction cor¬ 
responding to r G 7^ \ TZu- Since lo ■ {y!^ — yr) = 0 and p{y'j. — yr) = 0, we deduce that 
C(w) • [y’r-Vf) = 0. 

Let rv € 7^.2 for cr € A and rj = (^^(e) ~yr{e))- By Remark [T^ p{r]) = 0 and by definition 

C iv) = y'r ~ yra ■ Thus we have 


eScr 


C(w) • - yr„) = C(w) • C(?/) + p{^) ■ pih) =uj-r] = 0. 


This proves dTj). 
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m Assume that each connected component of Gu is strongly connected. Given oj € S’"*", 
we want to prove that there exist ai,, am such that 

LO — (wi , . . . , iOp, Q;i, . . . , CTm) S S . 

That is, such that coA = 0, where A is the stoichiometric matrix of {C, TZ) (see Subsection 17.ip . 
Let i = \R.u\- We order the set TZ in such a way that TZ \ TZu = ■ ■ ■■, fi}- Then, A can 

be written in block form as 

( ^ with Ai € A 2 and As € 

V ^2 ^mxie-I) J 

The columns of A 3 correspond to the reactions in TZi and are thus vectors of S by Definition 

ESI Hence, w [ „ ^ ~ choice of ai,..., am- It follows that 

\ J 


ujA 


0 'w’ {oj \; ■ ■ ■ ? ^p)A\ + (o^i, • • •, (y.Yfi^A2 — 0 


A 



1 

1 

) 



^■C{y'r. -Vrif 


^ ■ C,{y'rj-yrj) , 


= : V. 


(34) 


We can further reorder the species in U and the reactions in TZu, such that A 2 is a block 
diagonal matrix, where each block corresponds to one connected component of Gu ■ Due to this 
block structure, system ([5^ decomposes into subsystems given by the connected components 
of Gu- Therefore, it is enough to prove the existence of a for a strongly connected graph Gu- 
Hence we assume Gu is strongly connected. 

Let Cg be the incidence matrix of Gu and H as defined above Lemma [ 2 j We prove that 
Hv = 0. For fj € A, we have that & TZ and hence 


0 = UJ-{y~^- Vr,) = Ciy'r^e) “ Trie))- 

eG<T 


For a € A'\ A, strong connectedness and Remark EH imply that r((T) ^ 0. Thus, by definition 
of A, 

0 = ^ C {y'rie) - yr{e )) and hence 0 = ^ W • C (y(.(e) - yr{e)} - 

eGo- eGcr 


This shows that v € ker IT and hence by Lemma [2] we have u G im Cg. 

If Gu does not contain *, then A2 = Cg. Since v € im Cg we deduce that system (fMP has 
a solution (ai,... ,am)- If Gu contains *, then A 2 = Cg, with Cg being the matrix obtained 
from Cg by removing the last row, corresponding to the node *. We can then rewrite the 
system of equations (l3H) as 



( al^ 


/ai\ 

v = A\ 

) 

II 

VO/ 


(35) 


Since v G im Cg, there exists a vector b = ( 61 ,..., 6 ^+ 1 )* in ]^m+i Since 

the column sums of Cg are all zero by definition, (1,..., 1)* G ker Cg. Thus 


V = C^g{b-hm+l{l,---,lf)- 


Therefore, ai = bi — bm+i, for i = 1..., m is a solution to ([3^ . 


□ 
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9.5 Iterative elimination 


In this section we prove Proposition [22J We start with an auxiliary lemma. 

Lemma 3. Let lA he a linearly eliminable set such that Qu is connected and contains *. Let 
hi = T-LiU 'H 2 be a decomposition of hi such that U Let Qi{N) be the set of 

spanning trees rooted at N {*} in the sub-multidigraph of Qu- For any set of edges 

W in Gnu holds 

qu{x) ^ T^{T) = qni{x) ^ 7r(T). 

T&e{N),w<ZT Teei(v),iycT 

Proof. Let 02(*) be the set of spanning trees of Gh 2 rooted at *. Let r € 0(A^) be a spanning 
tree of Gu rooted at N. The path from any node N' G 'H 2 to N contains * by hypothesis. 
Therefore, r is the union of a spanning tree ti of Gni rooted at N and a spanning tree T 2 of 
Gh 2 rooted at *. Reciprocally, the union of any pair of spanning trees ti G 0i(A^), T 2 G 02(*) 
is a spanning tree of Gu rooted at iV. As subgraphs of Gu, and T 2 intersect at *. Hence 
7r(r) = 7r(Ti)7r(r2) and we obtain 

7r(r) = j Y 

rG0(Af),VrCT \riGei(V),VKCri / \T2Ge2(*) 

where we use that W is contained in Ghi- Using this computation with VU = 0 and N = *, 
we also have that 



qu{x) Y Y ^('^ 2 ) =qni{x) M Y ^('^ 2 ) 


kriG0i(*) 


i T2e02(*) 


kT2e02(*) 


Using these expressions the statement of the lemma follows. 


□ 


Note that we can write 

n(cj) = 7r(e) Y ^(7\e) 

7er{o-) 

for any edge e in a. The sum is over all spanning trees rooted at s(e) that contain the edges 
of (T\e. Lemma [3] guarantees that the computation of ^p{x) and the rate functions qu{x)Ii[a) 
of the reduced reaction network is independent of whether we consider Gu-, or Gni and Gn 2 
separately. Proposition [32] now follows because the sets of cycles of Gni and Gn 2 are disjoint. 


9.6 Post-translational modification networks 

Proof of Proposition |77[' PTM networks. The only edges in Gu whose associated reaction 
involves a substrate (that is, a species in IF^) are those with E as source or target. 

<^=) Let + E —... —>• Yj^ —>■ PI + 5*2, s > 0, be a path in (C, TZ). This path 
defines a cycle a in Gu with nodes E, ,Yj^ (it is a self-edge if s = 0). This cycle is 

such that r(cr) 7 ^ 0 (u belongs to a strongly connected component) and gives rise to a reaction 
Sn — 5 'j 2 in 112 . 

If Sn —>■ Si 2 G 7 ^ 2 ) then there is a cycle u G A defining it. The cycle a must contain 
a unique enzyme U G P. It follows that one of the edges in a corresponds to a reaction 
with reactant + E, and one of edges corresponds to a reaction with product Si^ + E. All 
other edges correspond to reactions between intermediates. The reactions corresponding to 
the edges in the cycle give the claimed path from 5*^ + PI to 5*2 + E. □ 
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